Plot list (for 3D plots scroll down from the legend; you can rotate and zoom the 3D plots with your mouse):
Volume Category Histograms | Where muons die (coffee ring plots) | Muon ranging out in iron yoke | Muon scatters in inflector | Distance distribution in Ring iron | Energy loss in volumes | Number of steps in volumes | Another muon in the iron | Stored muon | Coherent Betatron Oscillation
The goal of this document is to examine “escaped muons” that do not decay in the storage region and do not enter a calorimeter. Some of these muons may leave the storage ring altogether.
I generated art files of simulated events using mdc2 with different magnetic fields turned on or off (scenarios). Note: These files are still not the latest version of MDC2. Nathan is tuning up the injection gun. So right now the storage efficiency I see (0.5%) is still low. Should get better with the next round. We want to examine,
# Load necessary libraries
library(reticulate) # Access to python
library(stringr) # string functions
library(parallel) # Parallel processing (built-in to R)
library(dplyr) # data analysis
library(purrr) # Functional programming
library(readr) # Read formats
library(ggplot2) # ploting
library(rgl) # 3D plots
library(gridExtra) # Arranging plots
library(testthat) # testing
#
library(readGallery) # Read art Gallery files
# Put all readGallery::useDataProduct calls here
useDataProduct('std::vector<gm2truth::TrackingActionArtRecord>')
useDataProduct('std::vector<gm2ringsim::GeantTrackRecord>')
useDataProduct('std::vector<gm2truth::GhostDetectorArtRecord>')
# To restore most of the environment
trackingActionDF <- read_rds('trackingActionDF.rds')
trackingActionNDF <- read_rds('trackingActionNDF.rds')
volNameDF <- read_rds('volNameDF.rds')
ghCylDF <- read_rds('ghCylDF.rds')
ghWldDF <- read_rds('ghWldDF.rds')
Samples are located on the Fermilab dCache in the directory /pnfs/GM2/scratch/users/lyon/arr_20170321.
# Function to properly alter paths to accomodate XRootD
# e.g. convert /pnfs/BLA --> root://fndca1.fnal.gov/pnfs/fnal.gov/usr/BLA
xrootd_ify <- function(aPath) paste0('root://fndca1.fnal.gov/pnfs/fnal.gov/usr', str_replace(aPath, '/pnfs', ''))
Determine the locations of the data files
# Determine locations of our files
system("ifdh ls '/pnfs/GM2/scratch/users/lyon/arr_20170321/*/*.root' | grep ARR", intern=T) %>% sort() -> arrFiles
arrFiles %>% xrootify() -> arrXrdFiles
arrXrdFiles
[1] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243525_0/ARR_unified_everything_cyl.root"
[2] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243526_0/ARR_unified_noDipole_cyl.root"
[3] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243527_0/ARR_unified_noInflector_cyl.root"
[4] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243528_0/ARR_unified_noKickerQuads_cyl.root"
[5] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243530_0/ARR_unified_noKicker_cyl.root"
[6] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243532_0/ARR_unified_noQuads_cyl.root"
[7] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243540_0/ARR_unified_everything_cyl.root"
[8] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243542_0/ARR_unified_noDipole_cyl.root"
[9] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243544_0/ARR_unified_noInflector_cyl.root"
[10] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243547_0/ARR_unified_noKickerQuads_cyl.root"
[11] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243548_0/ARR_unified_noKicker_cyl.root"
[12] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243552_0/ARR_unified_noQuads_cyl.root"
write_rds(arrFiles, 'arrFiles.rds')
Pull out the scenaios (these will be in the same order as the file names)
scenarios <- str_match(arrFiles, 'arr_20170321/.+unified_(.+)_cyl')[,2]
scenarios %>% unique()
[1] "everything" "noDipole" "noInflector" "noKickerQuads" "noKicker" "noQuads"
Let’s look at the contents of these files (they’re all the same, so we’ll just do one) gm2 v7_04_00 will have product_sizes_dumper.
stringToRun <- str_interp(
"ssh lyon@gm2gpvm04.fnal.gov 'source /cvmfs/gm2.opensciencegrid.org/prod7/g-2/setup
setup gm2 v7_04_00 -q prof ;
product_sizes_dumper ${aFile} | grep ::'",
list(aFile = arrFiles[1] %>% str_replace('fnal.gov/usr/',''))
)
#
system(stringToRun)
1410762186 0.456 gm2truth::MagnetIronAndCryoArtRecords_artg4_magnetIronAndCryostat_mdc2arr.
660570363 0.214 gm2truth::TrackingActionArtRecords_artg4__mdc2arr.
637307490 0.206 gm2truth::TrajectoryArtRecords_artg4__mdc2arr.
240738334 0.078 gm2ringsim::GeantTrackRecords_artg4_KeepStepsAction_mdc2arr.
118035827 0.038 gm2truth::RingArtRecords_artg4_Ring_mdc2arr.
12078805 0.004 gm2truth::GhostDetectorArtRecords_artg4_GhostCylinderDetector_mdc2arr.
4482185 0.001 art::RNGsnapshots_randomsaver__mdc2arr.
4086898 0.001 gm2truth::RingTrackingPlaneArtRecords_artg4_RingTrackingPlanes_mdc2arr.
1373178 0.000 gm2truth::GhostDetectorArtRecords_artg4_GhostNearWorldDetector_mdc2arr.
1185772 0.000 gm2truth::PhaseSpaceArtRecord_artg4__mdc2arr.
36971 0.000 art::TriggerResults_TriggerResults__mdc2arr.
The tracking action data has birth and death data for every Geant track.
We’ll capture this information for the primary muon and its first generation daughters. See gm2dataproducts/mc/actions/track/TrackingActionArtRecord.hh.
We need a readGallery reader python class. Generate with,
readerClassSkel('gm2truth::TrackingActionArtRecord', writeFile = 'trackingActionReader.py')
Here is the reader,
read_file('trackingActionReader.py') %>% cat
from readGallery import GalleryReaderBase # Necessary for the base class
class TrackingActionArtRecordReader(GalleryReaderBase):
def __init__(self, inputTag):
GalleryReaderBase.__init__(self, inputTag)
self.names = ['fileEntry', 'eventEntry', 'trackType', 'trackID',
'parentTrackID', 'turn', 'volumeUID', 'status', 'p', 'e', 'x', 'y', 'z', 'px', 'py', 'pz']
def prepare(self, ROOT, ev):
GalleryReaderBase.prepare(self, ROOT, ev)
self.getValidHandle = ev.getValidHandle(ROOT.vector(ROOT.gm2truth.TrackingActionArtRecord))
def fill(self, ROOT, ev):
validHandle = self.getValidHandle(self.inputTag) # Get the valid handle for
# gm2truth::TrackingActionArtRecord
if not validHandle.empty(): # Does it have data?
p = validHandle.product() # Get the corresponding data product (maybe a vector)
# Fill from gm2truth::TrackingActionArtRecord
for e in p: # Loop over elements and fill
# Only accept the primary muon death and its first generation daughters birth
if (e.trackID == 1 and e.status == 1) or (e.parentTrackID == 1 and e.status == 0):
self.vals.append(
[ ev.fileEntry(), ev.eventEntry(), e.trackType, e.trackID,
e.parentTrackID, e.turn, e.volumeUID, e.status, e.p, e.e,
e.x, e.y, e.z, e.px, e.py, e.pz ])
return True
Make the reader class and object
trackingActionReaderC <- createReaderClass_from_file('trackingActionReader.py')$TrackingActionArtRecordReader
We have many files to process. Let’s do this reading in parallel.
We need to create a reader object per parallel job
trackingActionReaders <- map(1:length(arrXrdFiles), function(i) trackingActionReaderC(artInputTag("artg4")))
Load the data
# For some reason these jobs seem to run serially
jobs <- lapply(1:length(arrXrdFiles), function(i) getGalleryData(arrXrdFiles[i], trackingActionReaders[[i]]) %>%
mcparallel(name=i))
trackingActionReturns <- mccollect(jobs)
Let’s merge all of the output.
trackingActionDF <- map_df(1:length(trackingActionReaders),
function(i) galleryReader_df(trackingActionReaders[[i]]) %>%
mutate(scenario=scenarios[i],fileEntry=i))
nrow(trackingActionDF)
[1] 488953
trackingActionDF %>% head(100) # Let's not dump the entire data frame into the HTML page
# Write this out so we don't need to do the above again
write_rds(trackingActionDF, 'trackingActionDF.rds')
How many muons per file?
trackingActionDF %>% filter(trackID==1) %>% group_by(scenario, fileEntry) %>% tally()
We need to turn the volume IDs into volume names. The volume IDs change for each file. We can run a small art job to determine the volume ID to name tables.
runForVolIDString <- function(i) {
str_interp(
"PVS_CSVOUT=${csvout}_${i}_volNames.csv gm2 -c ${fclPath}/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 ${aFile}",
list( csvout=scenarios[i], i=i, fclPath=Sys.getenv("MRB_BUILDDIR"), aFile=arrXrdFiles[i])
)
}
runForVolIDStrings <- sapply(1:length(arrXrdFiles), runForVolIDString)
runForVolIDStrings
[1] "PVS_CSVOUT=everything_1_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243525_0/ARR_unified_everything_cyl.root"
[2] "PVS_CSVOUT=noDipole_2_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243526_0/ARR_unified_noDipole_cyl.root"
[3] "PVS_CSVOUT=noInflector_3_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243527_0/ARR_unified_noInflector_cyl.root"
[4] "PVS_CSVOUT=noKickerQuads_4_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243528_0/ARR_unified_noKickerQuads_cyl.root"
[5] "PVS_CSVOUT=noKicker_5_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243530_0/ARR_unified_noKicker_cyl.root"
[6] "PVS_CSVOUT=noQuads_6_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243532_0/ARR_unified_noQuads_cyl.root"
[7] "PVS_CSVOUT=everything_7_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243540_0/ARR_unified_everything_cyl.root"
[8] "PVS_CSVOUT=noDipole_8_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243542_0/ARR_unified_noDipole_cyl.root"
[9] "PVS_CSVOUT=noInflector_9_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243544_0/ARR_unified_noInflector_cyl.root"
[10] "PVS_CSVOUT=noKickerQuads_10_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243547_0/ARR_unified_noKickerQuads_cyl.root"
[11] "PVS_CSVOUT=noKicker_11_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243548_0/ARR_unified_noKicker_cyl.root"
[12] "PVS_CSVOUT=noQuads_12_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243552_0/ARR_unified_noQuads_cyl.root"
Let’s run the art jobs in parallel
jobs <- lapply(1:length(arrXrdFiles), function(i) mcparallel(system( runForVolIDStrings[i], intern=T ), name=i))
res <- mccollect(jobs)
and load the CSV files.
jobs <- lapply(1:length(arrFiles), function(i)
str_interp("${scn}_${i}_volNames.csv", list(scn=scenarios[i], i=i)) %>%
read_csv(col_names=c("volumeUID", "volName")) %>% mcparallel(name=i) )
#
volNameTables <- mccollect(jobs)
volNameDF <- map_df(1:length(volNameTables), function(i)
volNameTables[[i]] %>% mutate(scenario=scenarios[i],
fileEntry=i))
nrow(volNameDF)
[1] 1014
volNameDF %>% head(100)
write_rds(volNameDF, 'volNameDF.rds')
Now we merge with the tracking action data frame
trackingActionDF %>% inner_join(volNameDF) %>%
select(scenario, fileEntry, eventEntry, volName, everything()) -> trackingActionNDF
Joining, by = c("fileEntry", "volumeUID", "scenario")
nrow(trackingActionNDF)
[1] 488953
trackingActionNDF %>% head(100)
These muons decay in the storage region (not quite the same as “stored”),
trackingActionNDF %>% filter(trackID==1, volName=='ArcSection[00]') %>%
group_by(scenario, fileEntry) %>% tally()
Interesting! So, from this, only 0.575% decay in the storage region for the “all on” (or everything) sample. That seems rather low. But continuing.
write_rds(trackingActionNDF, 'trackingActionNDF.rds')
trackingActionNDF %>% filter(trackID == 1) %>% group_by(scenario, volName) %>% tally()
It may be nice to categorize the volume names by area (e.g. “Arc”, “bellows”, “Ring”, “Inflector”).
# Pull out the first word in camel case (see http://stackoverflow.com/questions/29916065/how-to-do-camelcase-split-in-python)
camelCaseSplit <- function(s) str_split(s, '(?<=[a-z])(?=[A-Z])|(?<=[A-Z])(?=[A-Z][a-z])')
# Test the above
expect_equal(camelCaseSplit("RingPoleTipUpper")[[1]], c("Ring", "Pole", "Tip", "Upper"))
expect_equal(camelCaseSplit("ringPoleTipUpper")[[1]], c("ring", "Pole", "Tip", "Upper"))
expect_equal(camelCaseSplit("world")[[1]], c("world"))
expect_equal(camelCaseSplit("ringPole")[[1]], c("ring", "Pole"))
# We really just want the first word
takeFirstCamel <- function(s) camelCaseSplit(s) %>% map_chr(function(s) s[1])
#
# Test
takeFirstCamel("RingPoleTipUpper") %>% expect_equal("Ring")
takeFirstCamel("ringPoleTipUpper") %>% expect_equal("ring")
takeFirstCamel("world") %>% expect_equal("world")
takeFirstCamel("ringPole") %>% expect_equal("ring")
takeFirstCamel(c("ringPoleTipUpper", "world")) %>% expect_equal(c("ring", "world"))
Don’t confuse RingTrackingPlane with the iron (also starts with Ring)
makeVolCategory <- function(s) ifelse( str_detect(s, 'RingTracking'),
"RngTrkPlane",
takeFirstCamel(s))
makeVolCategory(c("RingBottom", "RingTrackingPlane[8]", "world")) %>% expect_equal(c("Ring", "RngTrkPlane", "world"))
trackingActionNDF %>% mutate(volCategory = volName %>% makeVolCategory()) -> trackingActionNVDF
trackingActionNVDF %>% head(100)
Show count by category
trackingActionNVDF %>% filter(trackID == 1) %>% group_by(scenario, volCategory) %>% tally() -> volCatTally
volCatTally %>% mutate(perc = n/20000 * 100) -> volCatTally
volCatTally
trackingActionNVDF %>% filter(trackID ==1) %>% ggplot(aes(x=volCategory, group=scenario)) +
geom_bar()
ggplot(volCatTally, aes(x = volCategory, y=perc, fill=scenario)) +
geom_bar(stat="Identity", width=0.5, position="dodge")
Can’t really see non-Ring/world ones. Let’s split them out.
# Function to make this easy
plotit <- function(d) ggplot(d, aes(x = volCategory, y=perc, fill=scenario) ) +
geom_bar(stat="Identity", width=0.5, position="dodge") +
xlab("Volume category") + ylab("Percent")
volCatTally %>% filter(volCategory %in% c("Ring", "world")) %>% plotit -> p1
volCatTally %>% filter(! volCategory %in% c("Ring", "world")) %>% plotit -> p2
grid.arrange(p1, p2)
Let’s see where they die
# See https://rpubs.com/hadley/97970 for how to wrap a multipart ggplot2 component
plotCommon <- function () {
list(
facet_wrap(~scenario),
guides(col = guide_legend(override.aes = list(size=5, alpha=1))),
scale_color_discrete(name="Volume\nCategory"),
labs(x="z (mm)", y="x (mm)"),
theme_minimal()
)
}
trackingActionNVDF %>% filter(trackID == 1) %>%
ggplot( aes(x=z, y=x, color=volCategory)) +
geom_point(alpha=0.1) +
ggtitle('Where Muons Die') +
plotCommon()
Here it is without the ring losses, which dominate
trackingActionNVDF %>% filter(trackID == 1, volCategory != 'Ring') %>%
ggplot( aes(x=z, y=x, color=volCategory)) +
geom_point(alpha=0.5) +
ggtitle("Where Muons Die", subtitle = "Ring escapees excluded for clarity") +
plotCommon()
We have 3 GeV muons. Does it make sense that the vast majority of them stop in the iron? A table shows the Continuous Slow Down Approximation range of muons (CSDA). For a 3 GeV muon in iron, the CSDA is \(1.825 \times 10^{3}\) g/cm\(^2\). The density \(\rho\) is 7.874 g/cm\(^3\). So, where \(R\) is range,
\[R = \text{CSDA}/\rho\] Range is thus 2.32 m. This is bigger than the width of the iron, but remember that the muons come in at an oblique angle. We’ll have to prove this.
We have the Geant stepping information, so we should be able to figure out how far the muons go in the iron.
Here is the Gallery reader,
read_file('GeantTrackReader.py') %>% cat
from readGallery import GalleryReaderBase # Necessary for the base class
class GeantStepReader(GalleryReaderBase):
def __init__(self, inputTag, eventsToFill = []):
GalleryReaderBase.__init__(self, inputTag)
self.names = ['fileEntry', 'eventEntry', 'trackID', 'globalStepNum', 'totalEnergyDeposit',
'stepLength', 'volumeUID',
'x', 'y', 'z', 'px', 'py', 'pz']
self.eventsToFill = eventsToFill
print self.eventsToFill
def prepare(self, ROOT, ev):
GalleryReaderBase.prepare(self, ROOT, ev)
self.getValidHandle = ev.getValidHandle(ROOT.vector(ROOT.gm2ringsim.GeantTrackRecord))
def fill(self, ROOT, ev):
# Do we care about this event - if not then don't bother looking?
if self.eventsToFill and not ev.eventEntry() in self.eventsToFill:
return True
print 'Reading entry %s' % ev.eventEntry()
validHandle = self.getValidHandle(self.inputTag) # Get the valid handle for gm2ringsim::GeantTrackRecord
if not validHandle.empty(): # Does it have data?
p = validHandle.product() # Get the corresponding data product (maybe a vector)
# Fill from gm2ringsim::GeantStep
for e in p: # Loop over elements and fill
# Get the steps
for f in e.steps():
self.vals.append(
[ ev.fileEntry(), ev.eventEntry(), e.trackID(), f.globalStepNum(),
f.totalEnergyDeposit(), f.stepLength(), f.volumeUID(),
f.pos()[0], f.pos()[1], f.pos()[2],
f.p()[0], f.p()[1], f.p()[2] ])
# If we're on the last event, don't bother reading any more
if self.eventsToFill and ev.eventEntry() == self.eventsToFill[-1]:
return False
return True
geantStepReaderC <- createReaderClass_from_file('GeantTrackReader.py')$GeantStepReader
Let’s find an event where the muon did in the iron.
trackingActionNVDF %>% filter(trackID == 1, volName == 'RingYokeBottom') %>% head(100)
Event 4 looks good.
geantStepOne <- geantStepReaderC(artInputTag('artg4:KeepStepsAction'), tuple(as.integer(4)))
(4,)
getGalleryData(arrXrdFiles[1], geantStepOne)
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243525_0/ARR_unified_everything_cyl.root
Reading entry 4
Timings: Overall time = 5.234 s
Time per event: min=0.000 mean=0.206 max=1.027
gsOneDF <- galleryReader_df(geantStepOne)
Merge in volume names
gsOneDF %>% mutate(fileEntry = 1) %>%
inner_join(volNameDF) %>%
select(volName, everything()) -> gsOneDF
Joining, by = c("fileEntry", "volumeUID")
gsOneDF
Let’s make a 3D plot of this path
ring <- cylinder3d( center=rbind(c(0,-90, 0), c(0,90,0)),
radius=7112,
sides=20, closed = F)
Let’s do volume categories again.
gsOneDF %>% mutate(volCategory = volName %>% makeVolCategory %>% as.factor) -> gsOneDF
Let’s plot the path. Note that I have to make the legend separately since legend3d looks awful.
# lengend3d looks terrible - so do a regular legend
plot(1, type='n', axes=FALSE, ann=FALSE)
with(gsOneDF,
legend("top", legend=unique(volCategory), col=as.integer(unique(volCategory)), pch=19)
)
clear3d()
with(gsOneDF,{
plot3d(x=x, y=y, z=z, type='p', col = as.integer(volCategory), ylim=c(-400, 400))
})
plot3d(ring, add=T, alpha=0.2)
view3d(phi=90, theta=-90)
rglwidget()
How much energy loss and distance traveled in each volume category
gsOneDF %>% group_by(volCategory) %>% summarize(totalDist_meters=round(sum(stepLength)/1000, 2),
totalELoss_MeV=round(sum(totalEnergyDeposit), 1),
nSteps = n())
Ok - I believe this (note that the y scale is very small compared to the x,z scales - so the muon acutally travels quite far in x & z). The muon travels about 2.4m in the iron and ranges out. Physics (and Geant) work!
Let’s try another one
geantStepOneA <- geantStepReaderC(artInputTag('artg4:KeepStepsAction'), tuple(as.integer(10)))
(10,)
getGalleryData(arrXrdFiles[1], geantStepOneA)
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243525_0/ARR_unified_everything_cyl.root
Reading entry 10
Timings: Overall time = 2.566 s
Time per event: min=0.000 mean=0.090 max=0.993
gsOneADF <- galleryReader_df(geantStepOneA)
gsOneADF %>% mutate(fileEntry = 1) %>%
inner_join(volNameDF) %>%
mutate(volCategory = volName %>% makeVolCategory %>% as.factor)%>%
select(volName, volCategory, everything()) -> gsOneADF
Joining, by = c("fileEntry", "volumeUID")
gsOneADF
gsOneADF %>% group_by(volCategory) %>% summarize(totalDist_meters=round(sum(stepLength)/1000, 2),
totalELoss_MeV=round(sum(totalEnergyDeposit), 1),
nSteps = n()) -> gsOneADFSum
gsOneADFSum
Check sums
gsOneADFSum %>% summarize(sum(totalELoss_MeV))
gsOneADF %>% summarize(sum(totalEnergyDeposit))
Sanity check passes!
# lengend3d looks terrible - so do a regular legend
plot(1, type='n', axes=FALSE, ann=FALSE)
with(gsOneADF,
legend("top", legend=unique(volCategory), col=as.integer(unique(volCategory)), pch=19)
)
clear3d()
with(gsOneADF,{
plot3d(x=x, y=y, z=z, type='p', col = as.integer(volCategory), ylim=c(-400, 400), xlim=c(-10000, 10000),
zlim=c(-10000, 10000))
})
plot3d(ring, add=T, alpha=0.2)
view3d(phi=90, theta=-90)
rglwidget()
This one scatters in the inflector.
Can we show how far muons go in various materials?
Let’s collect data for 100 events.
We have a reader that can do n events.
read_file('GeantTrackReaderN.py') %>% cat
from readGallery import GalleryReaderBase # Necessary for the base class
class GeantStepReaderN(GalleryReaderBase):
def __init__(self, inputTag, startEvent=100, nEvents=10):
GalleryReaderBase.__init__(self, inputTag)
self.names = ['fileEntry', 'eventEntry', 'trackID', 'globalStepNum', 'globalTime', 'totalEnergyDeposit',
'stepLength', 'volumeUID',
'x', 'y', 'z', 'px', 'py', 'pz']
self.startEvent = startEvent
self.endEvent = startEvent + nEvents
def prepare(self, ROOT, ev):
GalleryReaderBase.prepare(self, ROOT, ev)
self.getValidHandle = ev.getValidHandle(ROOT.vector(ROOT.gm2ringsim.GeantTrackRecord))
def fill(self, ROOT, ev):
# Do we care about this event - if not then don't bother looking
if ev.eventEntry() < self.startEvent:
return True
if ev.eventEntry() >= self.endEvent:
return False
validHandle = self.getValidHandle(self.inputTag) # Get the valid handle for gm2ringsim::GeantTrackRecord
if not validHandle.empty(): # Does it have data?
p = validHandle.product() # Get the corresponding data product (maybe a vector)
# Fill from gm2ringsim::GeantStep
for e in p: # Loop over elements and fill
# Get the steps
for f in e.steps():
self.vals.append(
[ ev.fileEntry(), ev.eventEntry(), e.trackID(), f.globalStepNum(), f.globalTime(),
f.totalEnergyDeposit(), f.stepLength(), f.volumeUID(),
f.pos()[0], f.pos()[1], f.pos()[2],
f.p()[0], f.p()[1], f.p()[2] ])
return True
geantStepReaderNC <- createReaderClass_from_file('GeantTrackReaderN.py')$GeantStepReaderN
geantStepN <- geantStepReaderNC(artInputTag('artg4:KeepStepsAction'), 1000, 1000)
getGalleryData(arrXrdFiles[1], geantStepN)
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243525_0/ARR_unified_everything_cyl.root
Timings: Overall time = 51.662 s
Time per event: min=0.000 mean=0.024 max=0.452
gsNDF <- galleryReader_df(geantStepN)
gsNDF %>% mutate(fileEntry = 1,
p = sqrt(px*px + py*py + pz*pz)) %>%
inner_join(volNameDF) %>%
mutate(volCategory = volName %>% makeVolCategory %>% as.factor)%>%
select(volName, volCategory, p, everything()) -> gsNDF
Joining, by = c("fileEntry", "volumeUID")
nrow(gsNDF)
[1] 297501
gsNDF %>% head(100)
Let’s do a distribution of distance in the different volumes
gsNDF %>% group_by(eventEntry, volCategory) %>% summarize(
totalDist_meters=round(sum(stepLength)/1000, 2),
totalELoss_MeV=round(sum(totalEnergyDeposit), 1),
nSteps = n()) -> gsNDFSum
gsNDFSum
gsNDFSum %>% filter(volCategory == 'Ring') %>%
ggplot(aes(x=totalDist_meters)) + geom_histogram(bins=50) +
labs(x = 'Total distance in Ring material (m)', y='count', title='Muon in Ring Material')
Let’s look at this across the volume categories for energy loss
gsNDFSum %>% ggplot(aes(x=totalELoss_MeV)) + geom_histogram(bins=50) +
facet_wrap(~volCategory, scales = "free") +
labs(x = 'Total energy loss (MeV)', y='count', title='Energy loss in material',
subtitle = 'Note the different scales')
Let’s get an energy loss curve
gsNDF %>% ggplot(aes(x=p, y=totalEnergyDeposit)) + geom_point() + facet_wrap(~volCategory) +
labs(x="Muon Momentum (MeV/c)", y="Total Energy loss (MeV))")
Well, not sure what all this means. Let’s look for a particular event…
gsNDF %>% filter(eventEntry == 1300) %>% ggplot(aes(x=p, y=totalEnergyDeposit)) + geom_point() + facet_wrap(~volCategory) +
labs(x="Muon Momentum (MeV/c)", y="Total Energy loss (MeV))", title='Energy loss vs. Momentum for muons in event 1300' )
gsNDFSum %>% ggplot(aes(x=nSteps)) + geom_histogram(bins=50) +
facet_wrap(~volCategory, scales = "free") +
labs(x = 'Number of steps', y='count', title='Number of steps in materials',
subtitle = 'Note the different scales')
Let’s look at two events
gsNDFSum %>% filter(volCategory == 'Ring', totalDist_meters > 3)
Let’s look at event 1758
gsNDF %>% filter(eventEntry == 1758) -> gsNDF1758
gsNDF1758
gsNDF1758 %>% group_by(eventEntry, volCategory) %>% summarize(
totalDist_meters=round(sum(stepLength)/1000, 2),
totalELoss_MeV=round(sum(totalEnergyDeposit), 1),
nSteps = n()) -> gsNDF1758Sum
gsNDF1758Sum
Let’s plot it!
gsNDF1758 %>% mutate(volCategory = as.factor(as.character(volCategory))) -> gsNDF1758
# lengend3d looks terrible - so do a regular legend
plot(0, type='n', axes=FALSE, ann=FALSE)
with(gsNDF1758,
legend("bottom", legend=unique(volCategory), col=as.integer(unique(volCategory)), pch=19)
)
clear3d()
with(gsNDF1758,{
plot3d(x=x, y=y, z=z, type='p', col = as.integer(volCategory), ylim=c(-400, 400), xlim=c(-10000, 10000),
zlim=c(-10000, 10000))
})
plot3d(ring, add=T, alpha=0.2)
view3d(phi=90, theta=-90)
rglwidget()
Looks like this one started out in a strange place.
gsNDFSum %>% filter(volCategory == 'Arc', totalDist_meters > 870)
Let’s try event 1509
gsNDF %>% filter(eventEntry == 1509) -> gsNDF1509
gsNDF1509
gsNDF1509 %>% group_by(eventEntry, volCategory) %>% summarize(
totalDist_meters=round(sum(stepLength)/1000, 2),
totalELoss_MeV=round(sum(totalEnergyDeposit), 1),
nSteps = n()) -> gsNDF1509Sum
gsNDF1509Sum
Let’s plot it!
gsNDF1509 %>% mutate(volCategory = as.factor(as.character(volCategory))) -> gsNDF1509
# lengend3d looks terrible - so do a regular legend
plot(0, type='n', axes=FALSE, ann=FALSE)
with(gsNDF1509,
legend("bottom", legend=unique(volCategory), col=as.integer(unique(volCategory)), pch=19)
)
clear3d()
with(gsNDF1509,{
plot3d(x=x, y=y, z=z, type='p', col = as.integer(volCategory), ylim=c(-400, 400), xlim=c(-10000, 10000),
zlim=c(-10000, 10000))
})
plot3d(ring, add=T, alpha=0.2)
view3d(phi=90, theta=-90)
rglwidget()
Looks like a nice stored muon!! Can we see CBO?
gsNDF1509 %>% mutate(r = sqrt(x*x+z*z)) -> gsNDF1509
p1 <- gsNDF1509 %>% filter(volCategory %in% c('RngTrkPlane')) %>% ggplot(aes(x=globalTime, y=r)) + geom_line() +
labs(x='Global time (ns)', y='Radius (mm)', title='Radial CBO', subtitle="Measured by Ring Tracking Planes")
#
p2 <- gsNDF1509 %>% filter(volCategory %in% c('Arc','RngTrkPlane')) %>% ggplot(aes(x=globalTime, y=r)) +
geom_line() +
labs(x='Global time (ns)', y='Radius (mm)', subtitle="Measured by Arc hits and Ring Tracking Planes")
#
grid.arrange(p1, p2)
We can even see the kick!
Here we’ll plot the points too.
p1 <- gsNDF1509 %>% filter(volCategory %in% c('RngTrkPlane')) %>% ggplot(aes(x=globalTime, y=y)) + geom_line() +
geom_point() +
labs(x='Global time (ns)', y='y (mm)', title='Vertical CBO', subtitle="Measured by Ring Tracking Planes")
#
p2 <- gsNDF1509 %>% filter(volCategory %in% c('Arc','RngTrkPlane')) %>% ggplot(aes(x=globalTime, y=y)) +
geom_line() + geom_point() +
labs(x='Global time (ns)', y='y (mm)', subtitle="Measured by Arc hits and Ring Tracking Planes")
#
grid.arrange(p1, p2)
We want to examine muon escapees that leave the vacuum chamber. To measure this feature, we have a ghost cylinder at the outer wall of the vacuum chamber as well as above and below the storage region (before the muon would hit iron). That should catch every muon that doesn’t get stored and makes it out of the storage region. The data is gm2truth::GhostDetectorArtRecords_artg4_GhostCylinderDetector.
There is also a ghost cylinder just on the inside of the world cube to see muons that would actually leave the ring (and make it through the iron). The data is gm2truth::GhostDetectorArtRecords_artg4_GhostNearWorldDetector.
See gm2dataproducts/mc/ghostdetectors/GhostDetectorArtRecord.hh .
readerClassSkel('gm2truth::GhostDetectorArtRecord', writeFile = 'GhostDetectorReader.py')
read_file('GhostDetectorReader.py') %>% cat
from readGallery import GalleryReaderBase # Necessary for the base class
class GhostDetectorArtRecordReader(GalleryReaderBase):
def __init__(self, inputTag):
GalleryReaderBase.__init__(self, inputTag)
self.names = ['fileEntry', 'eventEntry', 'particleID',
'trackID', 'parentTrackID',
'x', 'y', 'z',
'px', 'py', 'pz',
'sx', 'sy', 'sz', 'time', 'volumeName']
def prepare(self, ROOT, ev):
GalleryReaderBase.prepare(self, ROOT, ev)
self.getValidHandle = ev.getValidHandle(ROOT.vector(ROOT.gm2truth.GhostDetectorArtRecord))
def fill(self, ROOT, ev):
validHandle = self.getValidHandle(self.inputTag) # Get the valid handle for gm2truth::GhostDetectorArtRecord
if not validHandle.empty(): # Does it have data?
p = validHandle.product() # Get the corresponding data product (maybe a vector)
# Fill from gm2truth::GhostDetectorArtRecord
for e in p: # Loop over elements and fill
self.vals.append(
[ ev.fileEntry(), ev.eventEntry(), e.particleID, e.trackID, e.parentTrackID,
e.position.x(), e.position.y(), e.position.z(),
e.momentum.x(), e.momentum.y(), e.momentum.z(),
e.spin.x(), e.spin.y(), e.spin.z(), e.time, e.volumeName ])
return True
ghostDetectorReaderC <- createReaderClass_from_file('GhostDetectorReader.py')$GhostDetectorArtRecordReader
Let’s read both ghost cylinders
ghCylReader <- ghostDetectorReaderC(artInputTag('artg4:GhostCylinderDetector'))
ghWldReader <- ghostDetectorReaderC(artInputTag('artg4:GhostNearWorldDetector'))
Load for one file
getGalleryData(arrXrdFiles, c(ghCylReader, ghWldReader))
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243525_0/ARR_unified_everything_cyl.root
Closing file, read 16740159 bytes in 100 transactions
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243526_0/ARR_unified_noDipole_cyl.root
Closing file, read 11886702 bytes in 66 transactions
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243527_0/ARR_unified_noInflector_cyl.root
Closing file, read 14680905 bytes in 66 transactions
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243528_0/ARR_unified_noKickerQuads_cyl.root
Closing file, read 16282559 bytes in 66 transactions
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243530_0/ARR_unified_noKicker_cyl.root
Closing file, read 16397132 bytes in 66 transactions
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243532_0/ARR_unified_noQuads_cyl.root
Closing file, read 16580366 bytes in 66 transactions
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243540_0/ARR_unified_everything_cyl.root
Closing file, read 16737803 bytes in 66 transactions
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243542_0/ARR_unified_noDipole_cyl.root
Closing file, read 11850233 bytes in 66 transactions
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243544_0/ARR_unified_noInflector_cyl.root
Closing file, read 14723528 bytes in 66 transactions
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243547_0/ARR_unified_noKickerQuads_cyl.root
Closing file, read 16377209 bytes in 66 transactions
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243548_0/ARR_unified_noKicker_cyl.root
Closing file, read 16397059 bytes in 66 transactions
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243552_0/ARR_unified_noQuads_cyl.root
Closing file, read 16814838 bytes in 66 transactions
Timings: Overall time = 252.037 s
Time per event: min=0.000 mean=0.000 max=7.044
Load the data into R
ghCylDF <- galleryReader_df(ghCylReader)
ghWldDF <- galleryReader_df(ghWldReader)
Look at what we got
nrow(ghCylDF) %>% format(big.mark=",")
[1] "2,374,697"
nrow(ghWldDF) %>% format(big.mark=",")
[1] "197,905"
ghCylDF %>% head(100)
ghWldDF %>% head(100)
write_rds(ghCylDF, 'ghCylDF.rds') # Write out just in case, since these are expensive to extract
write_rds(ghWldDF, 'ghWldDF.rds')
Fill in the scenario and pdg
pdgs <- c('e-'= 11, 'nu_e'= 12, 'mu-'= 13, 'nu_mu'= 14, 'tau-'= 15, 'nu_tau'= 16,
'e+'=-11, 'anti_nu_e'=-12, 'mu+'=-13, 'anti_nu_mu'=-14, 'tau+'=-15, 'anti_nu_tau'=-16,
'gam'=22, 'p'=2212, 'n'=2112, 'anti-p'=-2212, 'anti-n'=-2112)
pdgs
e- nu_e mu- nu_mu tau- nu_tau e+ anti_nu_e mu+ anti_nu_mu
11 12 13 14 15 16 -11 -12 -13 -14
tau+ anti_nu_tau gam p n anti-p anti-n
-15 -16 22 2212 2112 -2212 -2112
# We need to be a little careful, since the scenario names are repeated (we have two files per scenario).
# We can't make factor levels out of that.
# The file entries with 0-5 are easy. 6-11 are harder
addScenarioFactor <- function(f) factor(if_else(f >= 6, as.integer(f-6), f),
levels=0:(length(unique(scenarios))-1), labels=unique(scenarios))
#
mutateThis <- function(df) mutate(df,
scenario = addScenarioFactor(fileEntry),
pdg = factor(particleID, levels=pdgs, labels=names(pdgs)),
r = sqrt(x*x+z*z),
phi = atan2(z, x) %>% if_else(. > 0, ., 2*pi+ .) * (360.0/(2*pi)) )
#
ghCylDF %>% mutateThis -> ghCylDF
#
ghWldDF %>% mutateThis -> ghWldDF
The “ghost cylinder” virtual detector in the parallel world sees muon escapees that hit the vacuum chamber outer wall. How many hits per muon do we see?
Pull out just the muon
ghCylDF %>% filter(trackID == 1) -> ghCylMuDF
nrow(ghCylMuDF) %>% format(big.mark=',')
[1] "132,158"
To do this right, we need to make a “key” out of the fileEntry and eventEntry columns.
ghCylMuDF %>% mutate( flevKey = unite ... )
ghCylMuDF %>% group_by(fileEntry, eventEntry) %>% tally()
I’ve made 3D plots of all the scenarios (see below). Interesting to see what fraction get kicked up and down (are the ghost hits on the “ceiling or floor” or on the “wall”). We can tell by plotting the y coordinate…
ghCylDF %>%
ggplot(aes(y)) + geom_histogram(bins=50) + scale_y_log10() + facet_wrap(~scenario)
Let’s practice with just the everything scenario
ghCylMu_everything_DF <- ghCylMuDF %>% filter(scenario == 'everything')
nrow(ghCylMu_everything_DF) %>% format(big.mark=',')
[1] "23,691"
Let’s plot where these are
plotEscapees <- function(i) {
with( ghCylMuDF %>% filter(scenario == scenarios[i]),
plot3d(x=x, y=y, z=z, xlim=c(-8000, 8000), ylim=c(-120, 120), zlim=c(-8000, 8000), main=scenarios[i]))
plot3d(ring, add=T, alpha=0.2)
view3d(phi=90, theta=-90)
}
clear3d()
plotEscapees(1)
rglwidget()
They’re all lost up or down, not along the side. Let’s try more of the scenarios.
clear3d()
plotEscapees(2)
rglwidget()
No storage ring dipole field: As we expect, the muons don’t go very far.
clear3d()
plotEscapees(3)
rglwidget()
No inflector field: Seemingly lots of muons kicked up and down (need to check that). Many slam into the wall.
clear3d()
plotEscapees(4)
rglwidget()
clear3d()
plotEscapees(5)
rglwidget()
clear3d()
plotEscapees(6)
rglwidget()
Well, this isn’t too illuminating.